package scu.maqiang.peridynamics.interpolation;

import scu.maqiang.mesh.Tecplot;
import scu.maqiang.numeric.Algorithms;
import scu.maqiang.numeric.DirectDSolver;
import scu.maqiang.numeric.MVO;

import java.util.Arrays;

public class Interpolation1D {
    public static void main(String[] args) {
        //1. 输入参数
        double[] inputX = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
        double[] inputY = {2, 2, 2, 2, 2, 3.2, 3.2, 3.2, 3.2, 3.2};
        int inputSize = inputX.length;

        //2. 输出参数
        int outN = 200;
        double outDx = (inputX[inputSize - 1] - inputX[0]) / outN;
        double[] outX = MVO.linspace(0.5 * outDx, inputX[inputSize - 1] - 0.5 * outDx, outN);
        //System.out.println(Arrays.toString(outX));

        //3. 计算Lj_T
        double Eps = 1.0e-9;
        int p = 3;
        double[] Lj_T = new double[inputSize];
        double[] rho = new double[inputSize];
        for(int k = 0; k < outN; k++) {
            int kk = k;
            Arrays.setAll(rho, i -> 1.0/(Math.pow(Math.abs(outX[kk] - inputX[i]), 3) + Eps));
            double rhoSum = MVO.sum(rho);
            Arrays.setAll(Lj_T, i -> Lj_T[i] + rho[i] / rhoSum * outDx);
        }
        //System.out.println(Arrays.toString(Lj_T));

        double[] outY_N2 = new double[outN];
        double[] outY_N3 = new double[outN];
        double[] outY_N4 = new double[outN];
        double[] xi_kj = new double[inputSize];
        int[] order = new int[inputSize];
        int numberOfFamily = inputSize;
        double delta_k = 0.0;
        double[][] A_N2 = new double[3][3];
        double[][] A_N3 = new double[4][4];
        double[][] A_N4 = new double[5][5];
        double[] b_N2 = new double[3];
        double[] b_N3 = new double[4];
        double[] b_N4 = new double[5];
        double temp = 0.0;
        for(int k = 0; k < outN; k++) {
            int kk = k;
            MVO.fill(A_N2, 0.0);
            MVO.fill(A_N3, 0.0);
            MVO.fill(A_N4, 0.0);
            Arrays.setAll(xi_kj, i->Math.abs(outX[kk] - inputX[i]));
            Algorithms.selectionSort(xi_kj, order, true);
            delta_k = MVO.max(xi_kj);
            for(int jmem = 0; jmem < numberOfFamily; jmem++) {
                int j = order[jmem];
                double xi = inputX[j] - outX[k];
                double[] P_N2 = {1, xi, xi * xi};
                double[] weight_N2 = new double[3];
                temp = delta_k /(Math.abs(xi) + 1.0e-6);
                Arrays.fill(weight_N2, temp * temp);
                for(int ii = 0; ii < 3; ii++) {
                    for(int jj = 0; jj < 3; jj++) {
                        A_N2[ii][jj] += P_N2[ii] * P_N2[jj] * weight_N2[jj] * Lj_T[j];
                    }
                }

                double[] P_N3 = {1, xi, xi * xi, xi * xi * xi};
                double[] weight_N3 = new double[4];
                temp = delta_k /(Math.abs(xi) + 1.0e-6);
                Arrays.fill(weight_N3, temp * temp);
                for(int ii = 0; ii < 4; ii++) {
                    for(int jj = 0; jj < 4; jj++) {
                        A_N3[ii][jj] += P_N3[ii] * P_N3[jj] * weight_N3[jj] * Lj_T[j];
                    }
                }

                double[] P_N4 = {1, xi, xi * xi, xi * xi * xi, xi * xi * xi * xi};
                double[] weight_N4 = new double[5];
                temp = delta_k /(Math.abs(xi) + 1.0e-6);
                Arrays.fill(weight_N4, temp * temp);
                for(int ii = 0; ii < 5; ii++) {
                    for(int jj = 0; jj < 5; jj++) {
                        A_N4[ii][jj] += P_N4[ii] * P_N4[jj] * weight_N4[jj] * Lj_T[j];
                    }
                }
            }
            b_N2 = new double[]{1, 0, 0};
            b_N3 = new double[]{1, 0, 0, 0};
            b_N4 = new double[]{1, 0, 0, 0, 0};
            DirectDSolver solverN2 = new DirectDSolver(A_N2);
            double[] Acoef_N2 = solverN2.PGaussSolve(b_N2);
            DirectDSolver solverN3 = new DirectDSolver(A_N3);
            double[] Acoef_N3 = solverN3.PGaussSolve(b_N3);
            DirectDSolver solverN4 = new DirectDSolver(A_N4);
            double[] Acoef_N4 = solverN4.PGaussSolve(b_N4);
//            System.out.println(MVO.toString(Acoef_N2));
//            System.out.println(MVO.toString(Acoef_N3));
//            System.out.println(MVO.toString(Acoef_N4));
//            System.exit(1);
            for(int jmem = 0; jmem < numberOfFamily; jmem++) {
                int j = order[jmem];
                double xi = inputX[j] - outX[k];
                double yVal = inputY[j];

                double[] P_N2 = {1, xi, xi * xi};
                double[] weight_N2 = new double[3];
                temp = delta_k /(Math.abs(xi) + 1.0e-6);
                Arrays.fill(weight_N2, temp * temp);
                double[] tempVec = new double[P_N2.length];
                Arrays.setAll(tempVec, i -> Acoef_N2[i] * P_N2[i] * weight_N2[i]);
                double gfuncval_N2 = MVO.sum(tempVec);
                outY_N2[k] += yVal * gfuncval_N2 * Lj_T[j];

                double[] P_N3 = {1, xi, xi * xi, xi * xi * xi};
                double[] weight_N3 = new double[4];
                temp = delta_k /(Math.abs(xi) + 1.0e-6);
                Arrays.fill(weight_N3, temp * temp);
                tempVec = new double[P_N3.length];
                Arrays.setAll(tempVec, i -> Acoef_N3[i] * P_N3[i] * weight_N3[i]);
                double gfuncval_N3 = MVO.sum(tempVec);
                outY_N3[k] += yVal * gfuncval_N3 * Lj_T[j];

                double[] P_N4 = {1, xi, xi * xi, xi * xi * xi, xi * xi * xi * xi};
                double[] weight_N4 = new double[5];
                temp = delta_k /(Math.abs(xi) + 1.0e-6);
                Arrays.fill(weight_N4, temp * temp);
                tempVec = new double[P_N4.length];
                Arrays.setAll(tempVec, i -> Acoef_N4[i] * P_N4[i] * weight_N4[i]);
                double gfuncval_N4 = MVO.sum(tempVec);
                outY_N4[k] += yVal * gfuncval_N4 * Lj_T[j];
            }
        }
        System.out.println(Arrays.toString(outY_N2));
        System.out.println(Arrays.toString(outY_N3));
        System.out.println(Arrays.toString(outY_N4));
        Tecplot.LineXY("Interpolation1D.dat", outX, new double[][]{outY_N2, outY_N3, outY_N4});
    }

    public static double[] PD_Interpolation1D(double[] inputX, double[] inputY, double[] outX, int order) {
        int inputSize = inputX.length;

        //2. 输出参数
        int outN = outX.length;
        double outDx = (inputX[inputSize - 1] - inputX[0]) / outN;
        //3. 计算Lj_T
        double Eps = 1.0e-9;
        int p = 3;
        double[] Lj_T = new double[inputSize];
        double[] rho = new double[inputSize];
        for(int k = 0; k < outN; k++) {
            int kk = k;
            Arrays.setAll(rho, i -> 1.0/(Math.pow(Math.abs(outX[kk] - inputX[i]), 3) + Eps));
            double rhoSum = MVO.sum(rho);
            Arrays.setAll(Lj_T, i -> Lj_T[i] + rho[i] / rhoSum * outDx);
        }
        //System.out.println(Arrays.toString(Lj_T));

        double[] outY = new double[outN];
        double[] xi_kj = new double[inputSize];
        int[] distanceOrder = new int[inputSize];
        int numberOfFamily = inputSize;
        double delta_k = 0.0;
        double[][] A = new double[order + 1][order + 1];
        double[] b = new double[order + 1];
        b[0] = 1.0;
        double[] P = new double[order + 1];
        double[] weight = new double[order + 1];
        double temp = 0.0;
        for(int k = 0; k < outN; k++) {
            int kk = k;
            MVO.fill(A, 0.0);
            Arrays.setAll(xi_kj, i->Math.abs(outX[kk] - inputX[i]));
            Algorithms.selectionSort(xi_kj, distanceOrder, true);
            delta_k = MVO.max(xi_kj);
            for(int jmem = 0; jmem < numberOfFamily; jmem++) {
                int j = distanceOrder[jmem];
                double xi = inputX[j] - outX[k];
                Arrays.setAll(P, i -> Math.pow(xi, i * 1.0));
                temp = delta_k /(Math.abs(xi) + 1.0e-6);
                Arrays.fill(weight, temp * temp);
                for(int ii = 0; ii < 3; ii++) {
                    for(int jj = 0; jj < 3; jj++) {
                        A[ii][jj] += P[ii] * P[jj] * weight[jj] * Lj_T[j];
                    }
                }
            }
            Arrays.fill(b, 0.0);
            b[0] = 1.0;
            DirectDSolver solver = new DirectDSolver(A);
            double[] Acoef = solver.PGaussSolve(b);
            for(int jmem = 0; jmem < numberOfFamily; jmem++) {
                int j = distanceOrder[jmem];
                double xi = inputX[j] - outX[k];
                double yVal = inputY[j];
                double[] P_N2 = {1, xi, xi * xi};
                double[] weight_N2 = new double[3];
                temp = delta_k /(Math.abs(xi) + 1.0e-6);
                Arrays.fill(weight_N2, temp * temp);
                double[] tempVec = new double[P_N2.length];
                Arrays.setAll(tempVec, i -> Acoef[i] * P_N2[i] * weight_N2[i]);
                outY[k] += yVal * MVO.sum(tempVec) * Lj_T[j];
            }
        }
        return outY;
    }
}
